Take-Home Exercise 3

Author

Tan Wen Yang

Published

March 11, 2023

Modified

March 11, 2023

Loading in all the packages and data

NOTE: TAKE THE PROF’s MPSZ 2019 DATA FROM INCLASS-EX 09!!!! Rather than 2014

Compile to rds once the data wrangling part is done.

Geographically weighted methods (REgression / Random FoREST) and Ordinary Least Square Method (Conventional non-geo-weighted linear regression)

Note that corrplot cannot be loaded (Therefore has to be loaded from package)

pacman::p_load(sf, tidyverse, tmap, sfdep, onemapsgapi, httr, jsonlite, olsrr, ggpubr, GWmodel, dotenv, matrixStats)

Geospatial Data:

Singapore National Boundary and Master Plan 2019 subzone!

  • Singapore National Boundary is a polygon feature data showing the national boundary!

  • Master Plan 2019 subzone are information on URA 2019!

mpsz <- st_read(dsn="data/geospatial", "MPSZ-2019")
Reading layer `MPSZ-2019' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
national_boundary <- st_read(dsn="data/geospatial", "CostalOutline")
Reading layer `CostalOutline' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

Aspatial Data

  • Resale Flat prices

    • Training dataset from 1st January 2021 to 31st December 2022

    • Test dataset from 1 January to the last day of February 2023 resale prices For our

For our case, we will be using four-room flats for our prediction

resale_flat <- read_csv("data/aspatial/resale-flat-prices-from-jan-2017-onwards.csv")

Locational Factors:

For our assignment, we will need to look at location for

  • Proximity to CBD

  • Proximity to Eldercare

  • Proximity to Food court/Hawker

  • Proximity to MRT

  • Proximity to Park

  • Proximity to ‘Good’ primary schools

  • Proximity to Shopping Malls

  • Proximity to supermarket

  • Number of Kindergartens within 350m

  • Number of childcare centers within 350m

  • Number of primary schools within 1km

We start by sourcing some of this data

Geospatial Locational Source:

MRT/LRT and Bus Stop Locations

mrt_lrt <- st_read(dsn="data/geospatial", layer="RapidTransitSystemStation")
Reading layer `RapidTransitSystemStation' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 220 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 6068.209 ymin: 27478.44 xmax: 45377.5 ymax: 47913.58
Projected CRS: SVY21
bus_stop <- st_read(dsn="data/geospatial", layer="BusStop")
Reading layer `BusStop' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

There is something wrong, the MRT/LRT dataset is giving us polygon type rather than point. Despite the dataset saying it’s point.

OneMapSG Service API

We can extract the following according to their theme on OneMapSG API

  • Childcare (childcare)

  • Eldercare (eldercare)

  • Hawker Center (Queryname: hawkercentre)

  • Kindergarten (Queryname: Kindergartens)

  • Parks (Queryname: nationalparks)

Extra if we have time to think about:

  • Libraries (Queryname: libraries) [NUMBER WITHIN 350M]

  • Integrated Screening Programmes (queryname: moh_isp_clinics)

  • Tourist Attractions (Queryname: tourism) [NUMBER WITHIN 1KM]

Process of going through to create the shp file:

Courtesy of Megan’s work. (The following code chunk will be repeated to create the shp file for all the themes above).

#load_dot_env(file=".env")
#token <- Sys.getenv("TOKEN")
#themetibble <- get_theme(token, "themename")
#themetibble
#themesf <- st_as_sf(themetibble, coords=c("Lng", "Lat"), crs=4326)
#st_write(themesf, "themename.shp")

Load in all the shp data

childcare <- st_read(dsn="data/geospatial", layer="childcare")
Reading layer `childcare' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 1925 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
Geodetic CRS:  WGS 84
eldercare <- st_read(dsn="data/geospatial", layer="eldercare")
Reading layer `eldercare' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 133 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.7119 ymin: 1.271472 xmax: 103.9561 ymax: 1.439561
Geodetic CRS:  WGS 84
hawker_centre <- st_read(dsn="data/geospatial", layer="hawkercentre")
Reading layer `hawkercentre' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 125 features and 18 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6974 ymin: 1.272716 xmax: 103.9882 ymax: 1.449017
Geodetic CRS:  WGS 84
kindergarten <- st_read(dsn="data/geospatial", layer="kindergartens")
Reading layer `kindergartens' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 448 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6887 ymin: 1.247759 xmax: 103.9717 ymax: 1.455452
Geodetic CRS:  WGS 84
parks <- st_read(dsn="data/geospatial", layer="nationalparks")
Reading layer `nationalparks' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 421 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6929 ymin: 1.214491 xmax: 104.0538 ymax: 1.462094
Geodetic CRS:  WGS 84
libraries <- st_read(dsn="data/geospatial", layer="libraries")
Reading layer `libraries' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 31 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.7045 ymin: 1.263922 xmax: 103.9494 ymax: 1.448197
Geodetic CRS:  WGS 84
isp_clinics <- st_read(dsn="data/geospatial", layer="moh_isp_clinics")
Reading layer `moh_isp_clinics' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 378 features and 15 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6907 ymin: 1.26397 xmax: 103.9903 ymax: 1.456037
Geodetic CRS:  WGS 84
tourism <- st_read(dsn="data/geospatial", layer="tourism")
Reading layer `tourism' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 109 features and 15 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.684 ymin: 1.223 xmax: 103.974 ymax: 1.447
Geodetic CRS:  WGS 84

External Sources

  • Supermarkets (Source: Dataportal)

  • Primary School (Wikipedia, Reverse Geocoded using OneMap API)

  • Top Primary Schools (We will pick the top 10 based off SchLah’s dataset)

  • Malls (Taken from this dataset) [NOTE: SVY21]

supermarkets <- st_read("data/geospatial/supermarkets.kml")
Reading layer `SUPERMARKETS' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial\supermarkets.kml' 
  using driver `KML'
Simple feature collection with 526 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6258 ymin: 1.24715 xmax: 104.0036 ymax: 1.461526
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

Load Primary School Data with Lat/Long

The lat/long are taken from OneMapSG API

Top 10 schools as well, dervied from schlah

primary_sch <- read_csv("data/aspatial/primary_school_geo.csv")
top10_pri_sch <- read_csv("data/aspatial/top_10_primary_school_geo.csv")

Load Shopping Malls

shopping_mall <- read_csv("data/aspatial/shopping_mall.csv")

Data Wrangling (Geospatial)

List of task to do:

  • Convert some remaining data to Geospatial

  • Convert multipoint to point data (Removing the ‘z’) for supermarket

  • Convert all if not already in SVY21 into SVY21 (3414)

  • Remove all unnecessary columns

  • Check for null values

Convert all of our datasets from CSV to Geospatial (shp)

  • Primary Schools

  • Shopping Malls

Within primary school, the OneMapAPI is not able to properly return the correct coordinate as the name is ‘CATHOLIC HIGH SCHOOL (PRIMARY)’. For it to return the correct coordinate, the name has to be truncated without (PRIMARY).

primary_sch_sf <- st_as_sf(primary_sch, coords=c("LONG", "LAT"), crs=4326)
top_primary_sch_sf <- st_as_sf(top10_pri_sch, coords=c("LONG", "LAT"), crs=4326)
shopping_mall_sf <- st_as_sf(shopping_mall, coords=c("longitude", "latitude"), crs=4326)

Remove ‘Z-Dimension’ for Supermarket Data

supermarkets <- st_zm(supermarkets)

Check and change all EPSG for Geospatial data to ESPG 3414

childcare3414 <- st_transform(childcare, crs=3414)
eldercare3414 <- st_transform(eldercare, crs=3414)
hawker_centre3414 <- st_transform(hawker_centre, crs=3414)
kindergarten3414 <- st_transform(kindergarten, crs=3414)
parks3414 <- st_transform(parks, crs=3414)
libraries3414 <- st_transform(libraries, crs=3414)
isp_clinics3414 <- st_transform(isp_clinics, crs=3414)
tourism3414 <- st_transform(tourism, crs=3414)
primary_sch_sf_3414 <- st_transform(primary_sch_sf, crs=3414)
top_primary_sch_sf_3414 <- st_transform(top_primary_sch_sf, crs=3414)
shopping_mall_sf_3414 <- st_transform(shopping_mall_sf, crs=3414)

Check others:

st_crs(primary_sch_sf_3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(top_primary_sch_sf_3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(shopping_mall_sf_3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(supermarkets)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Looks like Supermarket is still in WGS84, let’s fix that

supermarkets3414 <- st_transform(supermarkets, crs=3414)

Check to make sure…

st_crs(supermarkets3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Remove all unnecessary columns

After glimpsing at the sf dataframe, we are going to remove all useless columns.

  • Information such as addresses/descriptions are not necessary, as all we need are the point data

  • We will keep the name for identification purposes

We will keep name and geometry

childcare3414 <- childcare3414 %>%
  select(c('NAME', 'geometry'))

Check result

childcare3414
Simple feature collection with 1925 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11810.03 ymin: 25596.33 xmax: 45404.24 ymax: 49300.88
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                             NAME                  geometry
1  APOLLO INTERNATIONAL PRESCHOOL PRIVATE LIMITED POINT (40985.94 33848.38)
2                    APPLE TREE PLAYHOUSE PTE LTD POINT (28308.65 45530.47)
3  Appleland Montessori Child Care Centre Pte Ltd POINT (17828.84 36607.36)
4                             APPLELAND PLAYHOUSE POINT (25579.73 29221.89)
5              APRICOT ACADEMY (LAGUNA) PTE. LTD. POINT (38981.02 32483.41)
6                                 Arise Preschool    POINT (21588.47 36307)
7        Artemis Preskool @ Tampines Pte Ltd (CC)  POINT (39239.78 37501.4)
8                    Artemis Preskool @ Woodleigh POINT (32389.52 35403.72)
9                     ARTS JUNIOR MONTESSORI  LLP POINT (25554.36 30090.82)
10                   Arts Kidz Pre-School Pte Ltd    POINT (28004.17 28442)

Same thing, name and geometry

eldercare3414 <- eldercare3414 %>%
  select(c('NAME', 'geometry'))

Check result

eldercare3414
Simple feature collection with 133 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                                           NAME
1                                  Yuhua Senior Activity Centre
2                                          THK SAC @ Kaki Bukit
3                                            THK SAC @ Boon Lay
4                        PEACE-Connect Senior Activity Centre@5
5                                        THK SAC @ Beo Crescent
6                                      Silver ACE @ Bukit Merah
7  Lions Befrienders Senior Activity Centre @ Tampines Blk 499C
8                    Care Corner Senior Activity Centre (WL569)
9           Fei Yue Senior Activity Centre (Bukit Batok Branch)
10       COMNET Senior Activity Centre @ 182 Rivervale Crescent
                    geometry
1  POINT (16614.08 36639.12)
2  POINT (38803.81 35098.77)
3  POINT (14481.92 36357.61)
4  POINT (31505.35 31853.52)
5  POINT (27218.35 30135.49)
6  POINT (27278.94 29350.17)
7  POINT (41665.14 37956.92)
8  POINT (23147.94 45761.17)
9  POINT (18820.58 36396.32)
10  POINT (36446.37 41376.9)
hawker_centre3414 <- hawker_centre3414 %>%
  select(c('NAME', 'geometry'))
hawker_centre3414
Simple feature collection with 125 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12874.19 ymin: 28355.97 xmax: 45241.4 ymax: 47850.43
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                   NAME                  geometry
1           Market Street Hawker Centre  POINT (29874.82 29650.7)
2          Marsiling Mall Hawker Centre POINT (22042.51 46139.03)
3          Margaret Drive Hawker Centre  POINT (24816.7 31094.91)
4       Fernvale Hawker Centre & Market  POINT (32867.9 41500.77)
5             One Punggol Hawker Centre POINT (35955.52 43336.13)
6          Bukit Canberra Hawker Centre POINT (26794.39 47850.43)
7                   Senja Hawker Centre POINT (19953.85 41008.06)
8                Buangkok Hawker Centre POINT (34575.37 40482.99)
9        Bukit Batok West Hawker Centre  POINT (17817.9 37478.45)
10 Telok Blangah Hawker Centre & Market     POINT (25191 28414.9)
isp_clinics3414 <- isp_clinics3414 %>%
  select(c('NAME', 'geometry'))
isp_clinics3414
Simple feature collection with 378 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12130.5 ymin: 27388.9 xmax: 45475.65 ymax: 48626.7
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                  NAME                  geometry
1             HEALTHWAY MEDICAL CLINIC POINT (35625.67 40809.19)
2             HEALTHWAY MEDICAL CLINIC  POINT (31090.4 33520.89)
3             HEALTHWAY MEDICAL CLINIC POINT (40903.43 39446.08)
4             HEALTHWAY MEDICAL CLINIC POINT (41384.47 37152.14)
5             HEALTHWAY MEDICAL CLINIC POINT (27371.14 45806.89)
6             HEALTHWAY MEDICAL CLINIC POINT (28014.68 45473.31)
7             HEALTHWAY MEDICAL CLINIC POINT (26462.45 47759.61)
8  HEALTHWISE MEDICAL CLINIC & SURGERY POINT (33725.38 32238.48)
9                    HL CLINIC PTE LTD  POINT (26451.91 28874.8)
10           HO MEDICAL CENTRE PTE LTD POINT (21103.55 32951.85)
kindergarten3414 <- kindergarten3414 %>%
  select(c('NAME', 'geometry'))
kindergarten3414
Simple feature collection with 448 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11909.7 ymin: 25596.33 xmax: 43395.47 ymax: 48562.06
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                                         NAME
1  PCF Sparkletots Preschool @ Cheng San-Seletar Blk 435 (KN)
2  PCF Sparkletots Preschool @ Cheng San-Seletar Blk 533 (KN)
3  PCF Sparkletots Preschool @ Cheng San-Seletar Blk 556 (DS)
4         PCF Sparkletots Preschool @ Chong Pang Blk 107 (KN)
5         PCF Sparkletots Preschool @ Chong Pang Blk 122 (KN)
6       PCF Sparkletots Preschool @ Chua Chu Kang Blk 10 (KN)
7           PCF Sparkletots Preschool @ Clementi Blk 330 (DS)
8              PCF Sparkletots Preschool @ Eunos Blk 616 (KN)
9           PCF Sparkletots Preschool @ Fengshan Blk 126 (DS)
10         PCF Sparkletots Preschool @ Fernvale Blk 416A (DS)
                    geometry
1  POINT (30325.45 38859.25)
2  POINT (30190.51 39574.18)
3   POINT (30705.05 39337.9)
4  POINT (27354.73 45992.92)
5  POINT (27755.87 46300.26)
6  POINT (19307.62 40271.08)
7  POINT (20706.39 32892.83)
8   POINT (37089.3 34892.34)
9  POINT (39752.13 34487.87)
10 POINT (33190.98 41392.37)
libraries3414 <- libraries3414 %>%
  select(c('NAME', 'geometry'))
libraries3414
Simple feature collection with 31 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 13665.24 ymin: 27383.57 xmax: 40922.89 ymax: 47759.75
Projected CRS: SVY21 / Singapore TM
First 10 features:
                           NAME                  geometry
1  Choa Chu Kang Public Library POINT (18175.64 40767.45)
2      Punggol Regional Library POINT (36014.52 43461.09)
3     Ang Mo Kio Public Library POINT (29364.62 39642.82)
4          Bedok Public Library POINT (38947.87 34357.48)
5         Bishan Public Library POINT (29727.41 36885.15)
6    Bukit Batok Public Library POINT (18646.83 36869.57)
7  Bukit Panjang Public Library POINT (20324.67 40211.98)
8        Central Public Library POINT (30327.79 31087.66)
9      Cheng San Public Library  POINT (34728.1 39397.71)
10      Clementi Public Library  POINT (20319.52 33038.5)

We only keep STN_NAM_DE and geometry

Note that the file is in polygon format, rather than point. Might have to deal with it.

https://stackoverflow.com/questions/52522872/r-sf-package-centroid-within-polygon

mrt_lrt <- mrt_lrt %>%
  select(c('STN_NAM_DE', 'geometry'))

Check

mrt_lrt
Simple feature collection with 220 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 6068.209 ymin: 27478.44 xmax: 45377.5 ymax: 47913.58
Projected CRS: SVY21
First 10 features:
                STN_NAM_DE                       geometry
1    ESPLANADE MRT STATION POLYGON ((30566.07 30621.21...
2   PAYA LEBAR MRT STATION POLYGON ((34495.6 33384.44,...
3  DHOBY GHAUT MRT STATION POLYGON ((29293.51 31312.53...
4       DAKOTA MRT STATION POLYGON ((34055.08 32290.62...
5     LAVENDER MRT STATION POLYGON ((31236.5 32085.76,...
6      RENJONG LRT STATION POLYGON ((34382.66 40949.64...
7        DOVER MRT STATION POLYGON ((21987.25 32576.91...
8      HOUGANG MRT STATION POLYGON ((34585.74 39132.47...
9      PHOENIX LRT STATION POLYGON ((19602.92 40048.64...
10    ALJUNIED MRT STATION POLYGON ((33422.35 33158.32...
parks3414 <- parks3414 %>%
  select(c('NAME', 'geometry'))
parks3414
Simple feature collection with 421 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12374.75 ymin: 21917.81 xmax: 52533.09 ymax: 49296.46
Projected CRS: SVY21 / Singapore TM
First 10 features:
                    NAME                  geometry
1      ENG KONG PLACE FC POINT (20647.02 35214.35)
2         JANGGUS GARDEN POINT (28409.23 48791.68)
3     JLN LIMAU MANIS PG POINT (40880.22 34107.22)
4         GARDEN VIEW PG POINT (31600.01 38049.84)
5       THOMSON GREEN PG POINT (27919.96 40074.49)
6           JLN RIANG PG POINT (31789.37 36584.45)
7   MEI HWAN CRESCENT PG POINT (31162.69 37212.14)
8          FULTON AVE PG   POINT (28026.6 38242.5)
9      MIMOSA TERRACE PG POINT (31016.83 40664.48)
10 JLN GENENG INTERIM PK  POINT (33295.7 37270.28)
primary_sch_sf_3414 <- primary_sch_sf_3414 %>%
  select(c('Primary school', 'geometry'))
primary_sch_sf
Simple feature collection with 181 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6878 ymin: 1.274958 xmax: 103.9628 ymax: 1.456608
Geodetic CRS:  WGS 84
# A tibble: 181 × 6
   `Primary school`         Type  Gender Locat…¹ TAG              geometry
 * <chr>                    <chr> <chr>  <chr>   <chr>         <POINT [°]>
 1 Admiralty Primary School Gove… Mixed  Woodla… ADMI…    (103.8 1.442635)
 2 Ahmad Ibrahim Primary S… Gove… Mixed  Yishun  AHMA… (103.8329 1.433153)
 3 Ai Tong School           Gove… Mixed  Bishan  AI T…  (103.833 1.360583)
 4 Alexandra Primary School Gove… Mixed  Bukit … ALEX… (103.8244 1.291334)
 5 Anchor Green Primary Sc… Gove… Mixed  Sengka… ANCH…  (103.8872 1.39037)
 6 Anderson Primary School  Gove… Mixed  Ang Mo… ANDE… (103.8414 1.384264)
 7 Ang Mo Kio Primary Scho… Gove… Mixed  Ang Mo… ANG … (103.8396 1.369322)
 8 Anglo-Chinese School (J… Gove… Boys   Central ANGL…   (103.841 1.30935)
 9 Anglo-Chinese School (P… Gove… Boys   Bukit … ANGL… (103.8356 1.318371)
10 Angsana Primary School   Gove… Mixed  Tampin… ANGS… (103.9518 1.348553)
# … with 171 more rows, and abbreviated variable name ¹​Location
supermarkets3414 <- supermarkets3414 %>%
  select(c('Name', 'geometry'))
supermarkets3414
Simple feature collection with 526 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4901.188 ymin: 25529.08 xmax: 46948.22 ymax: 49233.6
Projected CRS: SVY21 / Singapore TM
First 10 features:
     Name                  geometry
1   kml_1 POINT (35561.22 42685.17)
2   kml_2 POINT (32184.01 32947.46)
3   kml_3 POINT (33903.48 39480.46)
4   kml_4 POINT (37083.82 35017.47)
5   kml_5  POINT (41320.3 37283.82)
6   kml_6 POINT (41384.47 37152.14)
7   kml_7 POINT (30186.63 38602.77)
8   kml_8 POINT (28380.83 38842.16)
9   kml_9 POINT (34383.76 37311.19)
10 kml_10 POINT (29010.23 45755.51)
tourism3414 <- tourism3414 %>%
  select(c('NAME', 'geometry'))
tourism3414
Simple feature collection with 109 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11380.23 ymin: 22858.67 xmax: 43659.54 ymax: 47627.69
Projected CRS: SVY21 / Singapore TM
First 10 features:
                                                       NAME
1                                National Gallery Singapore
2                   Sultan Mosque (Masjid Sultan) Singapore
3           Sri Mariamman Temple: Hindu Temple in Singapore
4                              Armenian Church in Singapore
5                                         CHIJMES Singapore
6  St Andrewâ??s Cathedral- Singapore Architecture Landmark
7                                         Kreta Ayer Square
8                                  Albert Mall Trishaw Park
9                                     Chinatown Food Street
10                     Chinatown Heritage Centre, Singapore
                    geometry
1  POINT (30007.41 30267.17)
2  POINT (30877.14 31594.08)
3  POINT (29342.34 29382.57)
4   POINT (29818.66 30598.9)
5  POINT (30043.47 30820.05)
6  POINT (30113.58 30488.32)
7   POINT (29142.6 29293.89)
8  POINT (30244.32 31323.04)
9  POINT (29165.65 29443.75)
10 POINT (29227.71 29603.72)

Dealing with the anomaly - MRT/LRT

Despite it saying that’s point data, it gives polygon data, therefore we need to convert the polygon into point data. Unfortunately, there appears to be an issue that causes st_make_valid to not work, due to missing data within the polygon (There is an ‘na’).

Thus, we’ll take this chance to do some fine tuning and make use of OneMapSG to return us the correct coordinates. We will also remove all ‘depots’ as they are not considered MRT/LRTs for our dataset.

mrt_lrt_new_set <- read_csv("data/aspatial/mrt_lrt_lat_long.csv")

Convert into an SF Object

mrt_lrt_ns_sf <- st_as_sf(mrt_lrt_new_set, coords=c("LONG", "LAT"), crs=4326)
mrt_lrt_ns_sf <- mrt_lrt_ns_sf %>%
  select(c('STN_NAM_DE', 'geometry'))
mrt_lrt_ns_sf_3414 <- st_transform(mrt_lrt_ns_sf, crs=3414)

View the map to make sure the mapping is correct

tmap_mode("view")+
  tm_shape(bus_stop)+
  tm_dots(col="purple", size=0.05)

Looks like we have some bus stations in Malaysia, which makes sense, as there are some Singapore buses that travel into Malaysia.

These stops are:

  • LARKIN TERMINAL (Malaysia) [46239]

  • KOTARAYA II TER [46609]

  • JB SENTRAL [47701]

  • JOHOR BAHRU CHECKPOINT (46211, 46219)

We don’t want them in our calculation, so we’ll remove them.

bus_stop_to_remove <- c(46239, 46609, 47701, 46211, 46219)
bus_stop_cleaned <- filter(bus_stop, !(BUS_STOP_N %in% bus_stop_to_remove))

Check the map again:s

tm_shape(bus_stop_cleaned)+
  tm_dots(col="purple", size=0.05)

Convert Bus Stop from WGS84 to SVY21

bus_stop_cleaned_3414 <- st_transform(bus_stop_cleaned, crs = 3414)
st_crs(bus_stop_cleaned_3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Data Wrangling (Aspatial)

We will strictly be looking at Four Room Flats

  • Training Set: 1st January 2021 - 31st December 2022

  • Test Set: January - February 2023

resale_flat_sub <-  filter(resale_flat,flat_type == "4 ROOM") %>% 
              filter(month >= "2021-01" & month <= "2023-02")

Transforming the Data into usable data

Looking at the data, we will now need to combine both the block and street name to form the address so we can lookup for the postal code and retrieve the LAT/LONG

We will also need to calculate the remaining number of lease. Given that the remaining lease are in ‘YEAR’ + ‘MONTH’ format, we need to change it into a continuous number. We will change it to ‘MONTHS’.

rs_modified <- resale_flat_sub %>%
  mutate(resale_flat_sub, address = paste(block,street_name)) %>%
  mutate(resale_flat_sub, remaining_lease_yr = as.integer(str_sub(remaining_lease, 0, 2))) %>%
  mutate(resale_flat_sub, remaining_lease_mth = as.integer(str_sub(remaining_lease, 9, 11)))

Calculate the LEASE in Months and remove all intermediary columns

  • Convert all NA into 0

  • Sum up YEAR and MONTH and place them into a new column

  • Select only the columns we need.

rs_modified$remaining_lease_mth[is.na(rs_modified$remaining_lease_mth)] <- 0

rs_modified$remaining_lease_yr <- rs_modified$remaining_lease_yr * 12
rs_modified <- rs_modified %>% 
  mutate(remaining_lease_mths = rowSums(rs_modified[, c("remaining_lease_yr", "remaining_lease_mth")])) %>%
  select(month, town, address, block, street_name, flat_type, storey_range, floor_area_sqm, flat_model, 
         lease_commence_date, remaining_lease_mths, resale_price)

Glimpse at our data:

glimpse(rs_modified)

Retrieving the LAT/LONG Coordinates from OneMap.sg API

The following solution is courtesy of: Nor Aisyah

We first create a unique list of addresses

add_list <- sort(unique(rs_modified$address))

Now we will focus on getting the data needed.

In the following code chunk, the goal is to retrieve coordinates for a given search value using OneMap SG’s REST APIs. The code creates a dataframe to store the final retrieved coordinates and uses the httr package’s GET() function to make a GET request to https://developers.onemap.sg/commonapi/search.

The required variables to be included in the GET request are searchVal, returnGeom {Y/N}, and getAddrDetails {Y/N}.

The JSON response returned will contain multiple fields, but the code is only interested in the postal code and coordinates like Latitude & Longitude. A new dataframe is created to store each final set of coordinates retrieved during the loop, and the number of responses returned is checked to append to the main dataframe accordingly. Invalid addresses are also checked, and the response is appended to the main dataframe using rbind() function of base R package.

retrieve_coords <- function(add_list){
  
  # Create a data frame to store all retrieved coordinates
  postal_coords <- data.frame()
    
  for (i in add_list){
    #print(i)

    r <- GET('https://developers.onemap.sg/commonapi/search?',
           query=list(searchVal=i,
                     returnGeom='Y',
                     getAddrDetails='Y'))
    data <- fromJSON(rawToChar(r$content))
    found <- data$found
    res <- data$results
    
    # Create a new data frame for each address
    new_row <- data.frame()
    
    # If single result, append 
    if (found == 1){
      postal <- res$POSTAL 
      lat <- res$LATITUDE
      lng <- res$LONGITUDE
      new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
    }
    
    # If multiple results, drop NIL and append top 1
    else if (found > 1){
      # Remove those with NIL as postal
      res_sub <- res[res$POSTAL != "NIL", ]
      
      # Set as NA first if no Postal
      if (nrow(res_sub) == 0) {
          new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
      }
      
      else{
        top1 <- head(res_sub, n = 1)
        postal <- top1$POSTAL 
        lat <- top1$LATITUDE
        lng <- top1$LONGITUDE
        new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
      }
    }

    else {
      new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
    }
    
    # Add the row
    postal_coords <- rbind(postal_coords, new_row)
  }
  return(postal_coords)
}

We call the function to get the coordinates

coords <- retrieve_coords(add_list)

Inspect the results

coords[(is.na(coords$postal) | is.na(coords$latitude) | is.na(coords$longitude) | coords$postal=="NIL"), ]

From the results, we can see that it has missing postal code. However, the postal code is not as relevant to our analysis as latitude and longitude.

The missing respective postal code are:

  • 680215

  • 680216

rs_coords <- left_join(rs_modified, coords, by = c('address' = 'address'))

Write file to RDS

write_rds(rs_coords, "data/aspatial/rs_coords.rds")

Read the RDS file

resale_sub_flat <- read_rds("data/aspatial/rs_coords.rds")

Assign and transform CRS

Data from OneMapSG are in WGS84, as evident from the decimal. So we will create a SF object in that code before transforming into SVY21

resale_coords_sf <- st_as_sf(resale_sub_flat,
                    coords = c("longitude", 
                               "latitude"),
                    crs=4326) %>%
  st_transform(crs = 3414)
st_crs(resale_coords_sf)

Check if there are any NA outside of what we confirmed earlier

rows_with_na <- which(is.na(resale_coords_sf), arr.ind=TRUE)
if (nrow(rows_with_na) > 0) {
  message("The following rows have NA values:")
  print(rows_with_na)
} else {
  message("The dataframe does not contain any rows with NA values.")
}

Remove all unnecessary rows

We only need to keep the identifier, time, resale price, spatial point data, and the structural factors as listed below

  • Area of Unit

  • Floor Level

  • Remaining Lease

  • Age of Unit

trim_resale_flat <- resale_coords_sf %>%
  select(1:3, 6:8, 11:12)

Determining CBD

We need to consider the distance to CBD as well. We will take ‘Downtown Core’ as our reference point, which is located in the southwest of Singapore.

cbd_lat <- 1.287953
cbd_long <- 103.851784

cbd_sf <- data.frame(cbd_lat, cbd_long) %>%
  st_as_sf(coords = c("cbd_long", "cbd_lat"), crs=4326) %>%
  st_transform(crs=3414)

Proximity Distance Calculation

For our next step, we will need to integrate all of our geospatial data by calculating the proximity to the amenities.

The following function calculates proximity

  • rowMins from matrixStats package finds the shortest possible distance within the distance matrix.
proximity <- function(df1, df2, varname) {
  dist_matrix <- st_distance(df1, df2) %>%
    units::drop_units()
  df1[,varname] <- rowMins(dist_matrix)
  return(df1)
}

Implement the Proximity Calculations:

trim_resale_flat <- proximity(trim_resale_flat, cbd_sf, "PROX_CBD") %>%
  proximity(., eldercare3414, "PROX_ELDERCARE") %>%
  proximity(., hawker_centre3414, "PROX_HAWKER") %>%
  proximity(., isp_clinics3414, "PROX_ISP") %>%
  proximity(., mrt_lrt_ns_sf_3414, "PROX_MRT") %>%
  proximity(., parks3414, "PROX_PARKS") %>%
  proximity(., top_primary_sch_sf_3414, "PROX_TOP_PRI") %>%
  proximity(., shopping_mall_sf_3414, "PROX_SHOPPING") %>%
  proximity(., supermarkets3414, "PROX_SUPERMARKETS")

Facility Count within Radius

In addition to calculating the shortest distance between points, we’re also interested in finding out how many facilities are located within a specific radius. To accomplish this, we’ll use the st_distance() function to calculate the distance between the flats and the desired facilities. We’ll then use rowSums() to add up the observations and obtain the count of facilities within the desired radius. The resulting values will be added to the data frame as a new column.

num_radius <- function(df1, df2, varname, radius) {
  dist_matrix <- st_distance(df1, df2) %>%
    units::drop_units() %>%
    as.data.frame()
  df1[,varname] <- rowSums(dist_matrix <= radius)
  return(df1)
}

Implement the function:

trim_resale_flat <- 
  num_radius(trim_resale_flat, kindergarten3414, "NUM_KINDERGARTEN", 350) %>%
  num_radius(., childcare3414, "NUM_CHILDCARE", 350) %>%
  num_radius(., bus_stop_cleaned_3414, "NUM_BUS_STOPS", 350) %>%
  num_radius(., primary_sch_sf_3414, "NUM_PRI_SCHS", 1000) %>%
  num_radius(., tourism3414, "NUM_TOURIST_SITES", 1000)

Save the Data as RDS

Before we do that, we should trim the name, as when converting to RDS the column names are shortened due to limitations

trim_resale_flat <- trim_resale_flat %>%
  mutate() %>%
  rename("AREA_SQM" = "floor_area_sqm", 
         "LEASE_MTHS" = "remaining_lease_mths", 
         "PRICE"= "resale_price",
         "FLAT" = "flat_type", 
         "STOREY" = "storey_range")
st_write(trim_resale_flat, "data/geospatial/resale_flat_final.shp")

Reimport by reloading our dataset

resale_flat_sf <- st_read(dsn="data/geospatial", layer="resale_flat_final") 
Reading layer `resale_flat_final' from data source 
  `D:\Documents\IS415-GAA-WY\take-home-ex\take-home-ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 25503 features and 22 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11519.79 ymin: 28217.39 xmax: 42645.18 ymax: 48741.06
Projected CRS: SVY21 / Singapore TM

We’re done with gathering and creating our data! Now we can move onto our exploratory data analysis and building our prediction models. Note that there are still some data that need to be modified. Such as ordinal data like Storeys.

Speaking of which, we’ll handle them now.

Storey Data:

It’s a categorical data that has meaning when ordered. Higher/lower levels could potentially have an impact on the price of the HDB Flat. Therefore, we should convert this column of data into ordinal scale, one that the regression can learn from.